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ABSTRACT 

We consider the effectiveness of foreground cleaning in the recovery of Cosmic Microwave Background 
(CMB) polarization sourced by gravitational waves for tensor-to-scalar ratios in the range 0 < r < 0.1. Using 
the planned survey area, frequency bands, and sensitivity of the Cosmology Large Angular Scale Surveyor 
(CLASS), we simulate maps of Stokes Q and U parameters at 40, 90, 150, and 220 GHz, including realistic 
models of the CMB, diffuse Galactic thermal dust and synchrotron foregrounds, and Gaussian white noise. We 
use linear combinations (LCs) of the simulated multifrequency data to obtain maximum likelihood estimates 
of r, the relative scalar amplitude s, and LC coefficients. We find that for 10,000 simulations of a CLASS-like 
experiment using only measurements of the reionization peak (f < 23), there is a 95% C.L. upper limit of 
r < 0.017 in the case of no primordial gravitational waves. For simulations with r = 0.01, we recover at 
68% C.L. r = 0.012 +qqq*. The reionization peak corresponds to a fraction of the multipole moments probed 
by CLASS, and simulations including 30 < l < 100 further improve our upper limits to r < 0.008 at 95% 

C.L. (r = 0.010+OOQ4 for primordial gravitational waves with r = 0.01). In addition to decreasing the current 
upper bound on r by an order of magnitude, these foreground-cleaned low multipole data will achieve a cosmic 
variance limited measurement of the E-mode polarization’s reionization peak. 

Subject headings: cosmic background radiation - cosmological parameters - early universe - gravitational 
waves - inflation 


1. INTRODUCTION 

All astronomical data on cosmological scales conform to a 
six-parameter model of the Universe (ACDM) with dark mat¬ 
ter and a cosmological constant as the dominant components 
(Hinshaw et al. 2013; Planck Collaboration 2015c). The infla¬ 
tionary paradigm, which postulates a short period of exponen¬ 
tial expansion in the early Universe, accounts for features in 
ACDM, including the high degree of homogeneity and flat¬ 
ness and the scalar spectral index n s < 1 (e.g. Guth 1981; 
Starobinsky 1980; Kazanas 1980; Mukhanov & Chibisov 
1981; Einhorn & Sato 1981; Linde 1982; Albrecht & Stein- 
hardt 1982; Hinshaw et al. 2013; Bennett et al. 2013). One 
of inflation’s predictions is a super-horizon stochastic gravita¬ 
tional wave background that induces polarization in the CMB 
(Polnarev 1985). We can decompose the CMB’s polariza¬ 
tion field into E-modes with (-1)^ parity, and B-modes with 
( -l) M parity (Seljak & Zaldarriaga 1997; Kamionkowski 
et al. 1997; Hu & White 1997). An inflationary B-mode sig¬ 
nal can only come from primordial gravitational waves (tensor 
fluctuations of the metric), so a measurement of such a sig¬ 
nal would be strong evidence for an inflationary epoch. We 
quantify constraints on B-modes in terms of the ratio of ten¬ 
sor fluctuations to scalar (density) fluctuations, r, evaluated at 
k = 0.05 Mpc -1 . The current upper limit on tensor fluctua¬ 
tions (r < 0.09) comes from a combination of Planck and BI- 
CEP2 measurements (BICEP2/Keck et al. 2015; Planck Col- 
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laboration 2015e). 

The Cosmology Large Angular Scale Surveyor (CLASS) is 
a ground-based experiment that will observe 70% of the sky 
from Cerro Toco in the Atacama Desert at frequencies of 40, 
90, 150, and 220 GHz (Eimer et al. 2012; Essinger-Hileman 
et al. 2014; Rostem et al. 2014; Appel et al. 2014). First light 
for the experiment is imminent. CLASS’S 90 GHz band has a 
projected sensitivity of lOyuKarcmin, an improvement over 
the Planck 100 GHz map (118yuKarcmin at high-A Planck 
Collaboration 2015a), and will achieve the measurement sta¬ 
bility required to reach low multipoles using front end mod¬ 
ulation by a variable-delay polarization modulator (VPM) 
(Chuss et al. 2012). CLASS is currently the only planned 
sub-orbital mission exploring the combination of frequency 
and multipole ranges described above. CLASS will probe the 
reionization peak in the BB power spectrum (6 > 2°) along 
with other missions including PIPER (Lazear et al. 2014), 
QUIJOTE (Lopez-Caniego et al. 2014), LSPE (Aiola et al. 
2012), and GroundBird (Tajima et al. 2012). CLASS will 
also probe the recombination peak explored by BICEP2 (Og- 
burn et al. 2010), SPTpol (Austermann et al. 2012), ABS 
(Essinger-Hileman et al. 2010), ACTPol (Niemack et al. 
2010), POLARBEAR (Kermish et al. 2012), EBEX (Oxley 
et al. 2004), and SPIDER (Filippini et al. 2010). In contrast 
with other surveys that focus on higher multipoles or have 
different frequency ranges, CLASS probes a unique combina¬ 
tion of frequency and multipole space, as illustrated in Fig¬ 
ure 1. In addition to constraining inflation at all of the an¬ 
gular scales where B-modes are predominantly inflationary, 
we show that a noise-dominated BB spectrum Q oc has 
signal-to-noise that scales as 7 1-3 / 2 , as opposed to the well- 
known result t 111 found in the cosmic variance limit. In this 
era of initial measurements of B-mode polarization, we gain 
more information from large angular scale measurements than 
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Table 1 

Noise Contribution from Rescaled Foreground Templates 


Frequency v 
(GHz) 

-1/2 

Wy ' 

(jiK arcmin) 

a vt 

a s 

v -1/2 
a S W 40 
(ju K arcmin) 

ry v t 

v -1/2 
^220 

(ju K arcmin) 

40 

39 

1 

39 

0.022 

0.95 

90 

10 

0.103 

4.02 

0.095 

4.09 

150 

15 

0.032 

1.25 

0.306 

13.2 

220 

43 

0.018 

0.70 

1 

43 


Note. — For each band we list the expected 5-year polarized noise at each 
frequency (in“ 1/2 ). The columns for crfin“ 1/2 show the expected noise contri¬ 
bution in each band when we rescale the 40 and 220 GHz channels according 
to the synchrotron and dust template scalings used in simulations (Equation 2). 
The values of w ~ 1/2 are estimated in Essinger-Hileman et al. (2014). 

^ Values for the amplitude of synchrotron and dust templates at frequency v, orf, 
come from Equation 2 and assume fis = -3 and fin = 1.6. 

conventional wisdom would suggest, as we show in §A. 

In this paper we explore the recovery of B-modes in the 
presence of foregrounds, instrument noise, and a cut sky. This 
analysis will focus on a subset of the simulated CLASS data 
using only the multipoles € < 23, to distill the information 
available from measurements of the reionization peak alone. 
In §2 we describe simulations of the CMB with 0 < r < 0.1 
and foreground components in 40, 90, 150, and 220 GHz fre¬ 
quency bands using the CLASS sensitivity and sky coverage. 
In §3 we describe the maximum likelihood method for recov¬ 
ering the primordial B-mode signal in the presence of fore¬ 
grounds. We forecast constraints on r from this method ap¬ 
plied to the simulations in §4. 

2. MULTIFREQUENCY SIMULATIONS 

We use publicly available data to simulate the polarized 
Galactic synchrotron and dust foregrounds. Our synchrotron 
map templates are based on the WMAP9 K- band (23 GHz) 
polarization data (Bennett et al. 2013). Synchrotron polar¬ 
ization dominates the data in this band with negligible con¬ 
tribution from CMB polarization. Our dust polarization tem¬ 
plates are based on the Planck 353 GHz maps (Planck Col¬ 
laboration 2015a). In units of antenna temperature T A (v) = 
c 2 I v /2kBV 2 , with I y the specific intensity, we approximate the 
modified blackbody emission of dust as a power law with 
index p D = 1.6 (Planck Collaboration 2015f) and the syn¬ 
chrotron emission as a power law with fis = -3.0 (Ben¬ 
nett et al. 2013; Fuskeland et al. 2014). Using these power 
law indices we rescale the WMAP and Planck data to cre¬ 
ate templates for synchrotron and dust at 40 and 220 GHz, 
respectively. The minimum contamination from foregrounds 
is around 70 GHz (Bennett et al. 2013; Planck Collaboration 
2015b), which, along with the location of atmospheric water 
and oxygen lines, informs the CLASS experiment’s choice of 
band center frequencies (Eimer et al. 2012). 

We estimate CLASS’S sensitivity to tensor modes using 
random realizations of Gaussian CMB fluctuations with white 
Gaussian detector noise, all with the same model of fore¬ 
ground contamination. We use the fiducial ACDM param¬ 
eters from WMAP9 (Hinshaw et al. 2013) to simulate the 
scalar perturbation field, with a tensor contribution added in 
via Q = rCf nsor + sC/ alar . The contribution from gravita¬ 
tional lensing is included in the C^ calar term. We generate the 
theoretical spectra Q using CAMB, 5 smoothed in multipole 

5 http://camb.info 


space using a Gaussian window function with $fwhm = 15° 
to reduce aliasing from pixelization. We derive random Q 
and U maps from the theoretical Q values using the Healpix 
package synfast (Gorski et al. 2005) with pixels at resolu¬ 
tion #pi x ~ 7° (A^ide = 8), 6 including multipoles 2 < € <23. 
We then add dust and synchrotron polarization to each band, 
along with white Gaussian noise corresponding to the exper¬ 
iment sensitivity, given here in units of jjK arcmin (see Table 
1 ). 

Given CLASS’S wide scan strategy from its site in the Ata¬ 
cama Desert at latitude -23° (Essinger-Hileman et al. 2014), 
we include data with declination -73° < S < 27° so that 
/ sky ^0.7. We mask out the Galactic equator based on the 
WMAP P06 mask, defined in Page et al. (2007), resulting in a 
final / sky =* 0.5. 

We define the polarization maps P as a vector with the Q 
and U Stokes parameters as components. We use P s and P D 
as the synchrotron and dust templates at 40 GHz and 220 GHz 
respectively, and P CMB the simulated CMB. This allows us 
to generate the simulated multifrequency CLASS polarization 
data, 

P y = a v s P s + a v D P D + P CMB + N v (1) 


where N v contains the Gaussian white noise maps for Q and 
U, with the amplitude of the noise determined by the survey 
sensitivity given in Table 1 (Essinger-Hileman et al. 2014). 
Because we define the power law spectral indices Pi in terms 
of antenna temperature and our data are in units of thermo¬ 
dynamic temperature, we scale the foreground channels using 
the factor 


a y (h) = 


vp g(v)_ 
Vi) g(vd 


( 2 ) 


with g(v ) = dT/dT A = (e x - 1 ) 2 /(x 2 e x ) the conversion fac- 
tor from thermodynamic to antenna temperatures with x = 
hv/kTcMB = v/56.78 GHz. These a y give the relative strength 
of each foreground component i normalized by its strength 
at template frequencies v t . As described above these tem¬ 
plates are based on the WMAP K- band for synchrotron and 
Planck 353 GHz band for dust scaled to vs = 40 GHz and 
vd = 220 GHz, respectively. For these chosen frequen¬ 
cies and spectral indices, these strengths are = 0.103 
and = 0.095 at 90 GHz while for 150 GHz, we have 
o^ 50 = 0.032 and o^ 50 = 0.305, as shown in Table 1. We 
display sample maps of our simulated CLASS data in Figure 
2, along with the best-fit cleaned data and the residuals, to be 
described in §3. 

The spatial dependence of the spectral indices is based on 
the pre-flight Planck Sky Model (PSM) (Delabrouille et al. 
2013), which we display in Figure 3. As we demonstrate be¬ 
low, the PSM is consistent with current data. The Pi from this 
model come from intensity measurements of foreground com¬ 
ponents, and are not expected to be the same as the polarized 
power law indices. For this work, we rescale the PSM spec¬ 
tral indices to have mean values (J3s ) = -3.0 and (Pr>) =1.6. 
The amplitude of synchrotron spectral index fluctuations in 
the WMAP estimate derived from the K and Ka band data 


6 Healpix maps have resolutions denoted r = 0,1,2,.... Each of the 12 
lowest resolution pixels that characterize A^de = 1 divides into A^ide x Af s ide 
regions where /V s id e = 2 r , and the total number of pixels is N p \ x = 12A^ ide 
with characteristic pixel size 6 p \ x ~ 58.6° /Yide- The full documentation is at 
http://healpix.sf.net. 
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Angular scale A 6 ~ 180°/^ 

90 ° 10 ° 1 ° 0 . 1 ° 



£ range Foregrounds (ju K 2 ) 

Figure 1. We display schematically the regions of purview corresponding to balloon-borne and ground-based CMB polarization experiments. CLASS is unique 
in measuring both the reionization and recombination peaks while straddling the foreground minimum. As upper limits on r decrease, inflationary B-modes will 
dominate over lensing at increasing larger scales. The foreground model shown comes from Planck Collaboration (2015b) using measurements of 73% of the 
sky. We plot the frequency dependence of foregrounds in thermodynamic temperature units. 


40 GHz 90 GHz 



-6 6 -2 2 


150 GHz 220 GHz 



-4 /^ K 4 -12 12 


rfiihi dfljrx 

-0.5 /^ K 0.5 -0.5 0.5 

Figure 2. The top four plots show simulated CLASS observations of Q polar¬ 
ization assuming spatially varying spectral indices (as defined in Equation 2). 
The bottom two plots show (left) the output maximum likelihood CMB map 
with the simulation’s best-fit LC parameters as described in §3 and (right) the 
difference between the input CMB and the output maximum likelihood esti¬ 
mate. The difference map shows the effectiveness of our method for recon¬ 
structing the CMB, and the temperature range is representative of the noise 
in our map. Note that each of these maps uses a different color scale. We 
present the maps in Galactic coordinates (Mollweide projection) with gray 
pixels representing regions of the sky outside of the CLASS survey bound¬ 
ary or inside the Galactic mask. The black dashed lines on the 40 GHz map 
denote the declination limits of -73° and 27° in celestial coordinates. 

(Bennett et al. 2013; Fuskeland et al. 2014) is broadly con¬ 
sistent with the variation in the PSM’s /3s (Miville-Deschenes 


et al. 2008, model 4). The variation in the dust spectral index 
has variation A/3 D (Planck Collaboration 2015f) that is con¬ 
sistent with a constant spectral index (see §B) and hence with 
the small variation in the PSM (Finkbeiner et al. 1999, model 
8). Note that we only use the pre-flight PSM for spectral vari¬ 
ation, while we derive the 40 and 220 GHz polarization tem¬ 
plates from WMAP and Planck data respectively, as described 
above. However, the CLASS dataset itself will improve on 
this and be used in our final analysis. 

In this study we focus on the effects of foregrounds and 
a cut sky on the recoverability of primordial B-modes. To 
isolate these effects we do not include potential instrumental 
systematics. For measurements that use a VPM to suppress 
1 // noise, Miller et al. (2015) find that B-mode power from 
systematic effects can be suppressed to levels below the pri¬ 
mordial level for r = 0.01, which suggests such effects are 
negligible in the context of this work. 

3. PIXELIZED LIKELIHOOD ANALYSIS 

We base our analysis on methods developed in Efstathiou 
(2006) and Katayama & Komatsu (2011), and extend the 
Gaussian likelihood formalism to an arbitrary number of fre¬ 
quency channels, each with its own linear fit coefficient. We 
define the likelihood 

exp [-\x(a v ) T C~ l (r, s, a v )x(a v )] 

£ oc- - =-, (3) 

ydet[C p (r, s, a v )] 

where 

n n 

x = X a v pv > X Uv = 1 H) 

V V 

is the CMB map cleaned using n frequency channels, with 
£ v a v = 1 preserving the CMB’s frequency spectrum. We 
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Figure 3. The top row of foreground spectral indices from the PSM (De- 
labrouille et al. 2013) are used in our simulations. We model the variation 
in spectral index by breaking the sky into 14 regions (middle row, labelled 
“Mean”) with distinct LC fit coefficients. The fit in the middle row are the 
means of the PSM spectral indices within each region. This modelling trades 
off accuracy at small scales for computational efficiency. 


achieve this constraint defining ago as 1 - <240 - a 150 - < 2220 - We 
define the pixel-space covariance matrix, 

C p = rC* nsor + sC^ r + IpO 2 . (5) 

Here I p is the identity matrix in pixel space, and cr 2 = 
Yjv( a vO~v) 2 is the survey noise. The C p ’s are the theoretical 
signal covariance matrices, which themselves come from the 
power spectra, C BB and C EE , assuming uncorrelated E- and B- 
modes, as detailed in Appendix A of Katayama & Komatsu 
(2011). We vary s, the amplitude of scalar fluctuations nor¬ 
malized by a fiducial value, to mitigate spurious correlations 
between the CMB and foregrounds in the maximum likeli¬ 
hood fitting as in §5.1 of Katayama & Komatsu (2011). We 
use uniform unbounded priors on all parameters except for 
the tensor-to-scalar ratio and the 40 GHz and 220 GHz coef¬ 
ficients. We impose the prior r > 0 to avoid singular covari¬ 
ance matrices in Equation 5, and we impose the prior that 
the 40 and 220 GHz coefficients are negative to avoid nu¬ 
merical errors (which does not affect our overall results, as 
seen in Figure 4). The product of our likelihood and priors is 
proportional to the Bayesian posterior probability distribution 
for our parameters. We emphasize that this method differs 
from the internal linear combination method described in Ef- 
stathiou et al. (2009) due to our simultaneous variation of cos¬ 
mological and foreground LC parameters and also due to the 
presence of the C~ l weighting for the solution of LC coeffi¬ 
cients. When we maximize the likelihood (Equation 3) using 
the covariance matrix (Equation 5), the maximum likelihood 
LC coefficients, {a v }, we create a map, x , that best approx¬ 
imates a realization of the CMB with Gaussian noise. The 
modified foreground removal method is effective because the 
pixel-pixel correlations, C p , drive the fit so that the mean of 
the map, ( x ), is zero, with variance corresponding to CMB 
fluctuations and instrumental white noise. 


We obtain an estimate of the CMB, x , by finding the pa¬ 
rameters r, s, < 240 , < 290 , < 2 i 50 > and <2220 that maximize the like¬ 
lihood £ with the priors r > 0 and < 240 , (2220 < 0- We find 
that although a Monte Carlo Markov Chain (MCMC) esti¬ 
mation gives posterior distributions and suitable estimators 
for these parameters, using the scipy minimization routine 
fmin_l_bfgs_b, an implementation of the limited-memory 
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Byrd 
et al. 1995), gives estimates of the maximum likelihood val¬ 
ues for parameters {?ml, < 2 v,ml} at a fraction of the time 

required to compute the full MCMC posterior distributions. 
We plan to use MCMC methods for the analysis of the actual 
multi-frequency CLASS observed maps, but cosmic variance 
can cause different posterior distributions for the same under¬ 
lying value of r. For our ~ 10 4 simulations with random real¬ 
izations of CMB and noise, the distribution of maximum like¬ 
lihood values for r is comparable to the width of the posterior 
distribution of a single Monte Carlo chain. For forecasting 
purposes, we therefore run large numbers of simulations and 
study the posterior distribution of maximum likelihood values 
of the model parameters. 

We also model spatially varying spectral behavior in the 
foregrounds. In practice, we split the map into fourteen sep¬ 
arate foreground regions defined by the intersection of un¬ 
masked dust and synchrotron regions shown in Figure 3, 
which results in a 44-parameter fit (four a v coefficients with 
one constraint x 14 regions, s, and r) for varying spectral be¬ 
havior. We define the regions in the middle row of Figure 
3 using the PSM components in the top row. To construct 
these regions, we smooth the PSM pi maps, find all local min¬ 
ima and maxima, and choose contiguous pixels with roughly 
constant spectral index values. We choose thresholds such 
that the regions are small enough to make a constant index a 
good approximation, and large enough that we maintain com¬ 
putational efficiency that follows from the smaller number of 
parameters. Defining these regions according to existing mea¬ 
surements of spectral variation is the only part of our analysis 
that requires data external to CLASS. When the CLASS data 
are available, they will be the best probe of this large-scale 
position dependence of the foreground spectral indices. 

We discuss the distribution of maximum likelihood r values 
in the following section. For typical simulations with input 
r = 0 . 01 , the other parameters for a single representative re¬ 
gion have sample mean and standard deviation s = 0.97+0.13, 
<240 = -0.1 ± 0.01, ago - 0.9 ± 0.1, (2150 = 0.3 ± 0.1, and 
<2220 = -0.16 ± 0.03. Both (240 and (2220 are negative because 
data at these frequencies serve to remove foregrounds, while 
the 90 and 150 GHz coefficients contribute positively, with 
most of the weight coming from ago. We display the distribu¬ 
tion of ML values for these parameters and r in Figure 4. 

4. RECOVERY OF r 

If the uncertainty were only dependent on instrument noise 
and not a function of the cosmological value of r, we would 
choose a threshold based on this uncertainty, such as 3<x r , and 
claim this as a detection limit. However, the prediction prob¬ 
lem is more complicated, as cr r depends on r as well as dust 
and synchrotron emission in the Galaxy through chance cross¬ 
correlations. We focus on the effect of varying r on the detec¬ 
tion limit, which we will not know until real measurements 
are made. We find the distribution of recovered maximum 
likelihood values tml for a given input r. 

Using the simulations described in §2 and the maximum 
likelihood method described in §3, we generate realizations of 
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r ML ^ML ^40 a 90 a l50 a 220 


Figure 4. Distribution of the output maximum likelihood parameters for 10 4 simulations using the sensitivities from Table 1 with an input r = 0.01. The 
simulation features spatially varying spectral indices and the distribution of LC coefficients a v are displayed for one of the fourteen regions in Figure 3. This 
demonstrates that the maximum likelihood solutions for hvil and sml are uncorrelated with each other and the LC coefficients. 


CLASS-like data and infer best-fit parameters from them. The 
maps all have scalar fluctuations with a fixed normalization 
(s = 1) and the same foreground simulation from §2, while 
the tensor modes have varying strength 0 < r < 0.1 uniformly 
sampled with one r every Sr = 5 x 10 -6 . We display a subset 
of these simulations in Figure 5, marginalized over s and the 
LC coefficients. 

Actual CLASS data will provide a maximum likelihood es¬ 
timate tml- Using simulated data similar to those from Fig¬ 
ure 5, we will estimate the distribution of r using the slice 
p(A r ML)- Without knowledge of t M l, however, we need to 
marginalize over the possible values of tml that result from 
cosmic variance given a “true” r. 

We estimate the constraining power of CLASS prior to a 
specific measurement tml- Our goal is to calculate a proba¬ 
bility distribution p(r\r tme ) that indicates our best estimate of 
r given some true value, r tme . From simulations, we calculate 
the probability of a specific value of r given a measured r M l, 
p(r\r ml), and the probability of obtaining a given maximum 


likelihood result r ML given the true value of r, p(r M Lktme)- 
Using this information, we can marginalize over the nuisance 
parameter, r ML , i.e. 

POktrue) = J' P(r\rML)PiruhVtme) d' ML- (6) 

The quantity p(r\r tme ) represents the average of a CLASS-like 
experiment’s posterior distribution for r over many realiza¬ 
tions, given that the true amplitude of tensor fluctuations is 
r tme . This statistic is motivated by noting a single CLASS-like 
experiment will produce a probability distribution p(r|rML)- 
The random variable r M L is a draw from the distribution 
p(r MLktraeX and is not of inherent interest to this discussion. 
The p(rMLktme) function acts as a kernel for a weighted av¬ 
erage, and is used here to create a weighted average of each 
p(?\ r ML), which results in the distribution p(r\r tme ). We illus¬ 
trate this process in Figure 5 in which the histogram across the 
top gives a single realization of p(r|r M L) and the histogram on 
the right shows /?(r M Lktme) for a specific value r = r true for the 
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Figure 5. Each point in the bottom-left plot represents a unique realization 
of the simulations described in §2. Given an input tensor-to-scalar ratio r 
we compute the distribution of recovered maximum likelihood values hvil- 
We also display sample distributions of hvil given an r trU e = 0.01 ± 0.0005 
(green histogram) in the right plot, and of r given an hvil = 0.01+0.0005 (red 
histogram) in the top plot. To obtain our posterior distribution of r given a 
value r tme , we integrate the product of these distributions over hvil (Equation 
6). 

input model. In this formalism, we can obtain a distribution 
of expected estimates of r given r true . 

We also estimate the constraining power of CLASS’S high- 
multipole data, with 0 pix ~ 27.5' (N s ^ e = 64, £ < 100). 
The model used here has a fixed scalar fluctuation amplitude 
s = l and spatially constant spectral indices. These simplifi¬ 
cations are justified because foreground power spectra decay 
approximately as Q oc Z’ -2 - 5 , and spurious correlations be¬ 
tween the CMB and foregrounds are negligible at these scales 
(Page et al. 2007; Planck Collaboration 2014). We use an in¬ 
ternal linear combination (ILC) method similar to that used 
in the WMAP analysis (Hinshaw et al. 2007) to clean fore¬ 
grounds, and then use PolSpice (Chon et al. 2004) to extract 
the B-mode power spectrum for multipoles 30 < l < 100. 
By simulating maps with varying r as in the low-/’ simula¬ 
tions, we find that we can describe the likelihood of r from 
the high-/’ distribution with a Gaussian probability distribu¬ 
tion p(r\r tme ) oc £-( r - r true) 2 /2o- 2 ^ _ q.005. 

This analysis omits multipoles 24 < l < 29, due to con¬ 
straints from both the low -l and high-/’ software. For the 
pixel-based approach at low-G the Healpix documentation 
recommends using £ max = - 1 = 23 to avoid oversam¬ 

pling each pixel, while PolSpice introduces mode-coupling 
effects in the power spectrum for l < 30. 

With the distribution for r tme = 0, we obtain the range of 
potential false positive results, which suggests a definition of 
our upper limit as the value of r that contains 95% of the area 
under the curve from r = 0, yielding r < 0.017 for low-/’ data 
alone, and r < 0.008 using the high-/ 1 simulations as well. 
Similarly, we estimate the expected range of recovered r for 
rtme = 0.01, and find that the 68% confidence interval, de¬ 
fined as the percentiles (0.16,0.5,0.84), gives r = 0.01 
(r = 0.010+OOQ4 including high-/’). These distributions are dis¬ 
played in Figure 6. Using a simple likelihood ratio test based 
on these curves with Jif(r) = p(r\r tmQ = 0.01), we find a like¬ 
lihood ratio for the low-/ 1 data with J2f(0)/Jjf(0.01) = 0.52, 
corresponding to a PTE of 0.13, or a significance of 1.2<x. In- 



r 


Figure 6. We estimate the distributions for r from Equation 6, and use a 
Gaussian kernel density estimator to obtain a continuous curve, which re¬ 
turns a normalized probability density by default. The dashed lines are the 
distributions given only the low-^ data (2 < t < 23), while the solid curves 
are a product of the low-^ distributions and a Gaussian estimate of the high- 
/ (30 < t < 100) distribution with mean \i - r tr ue and standard devia¬ 
tion cr = 0.005. The solid curve for r tme = 0 gives a 95% C.L. upper 
limit r < 0.008, while for the other solid curves the 95% C.L. limits are 
0.006 < r < 0.014 and 0.015 < r < 0.024. 

eluding high-/’ data gives a projected Jjf(0)/Jjf(0.01) = 0.07, 
corresponding to a PTE of 0.01, or a significance of 2.3cr. 

5. CONCLUSIONS 

In this work, we have simulated CMB polarization maps on 
a cut sky contaminated by diffuse foreground emission and 
Gaussian white noise. We recovered the input CMB by solv¬ 
ing for maximum-likelihood maps using a pixel-based maxi¬ 
mum likelihood method tailored to measure the reionization 
peak (Z’ < 23). We have found the following. 

• We found unbiased estimates of the amplitudes of 
scalar and tensor fluctuations with a four-frequency 
maximum likelihood analysis accounting for CLASS - 
like noise, a cut sky, and foregrounds with spatially 
varying spectral indices. 

• We have found that, by using only polarization the 
reionization peak (2 < t < 23), data from a five-year 
CLASS-like experiment would produce a 95% upper 
limit of r < 0.017, which improves to r < 0.008 when 
we include high-/’ data (30 < € < 100). 

• Our pixel-based low-multipole maximum likelihood 
solution is straightforward to run and implement, and is 
stable to foreground complications, i.e. variable spec¬ 
tral indices and polarization fractions. 

Although CLASS’S primary purpose is constraining the 
amplitude of tensor fluctuations, a foreground-cleaned mea¬ 
surement of CMB polarization on the largest scales can help 
constrain several other cosmological parameters. Our cos¬ 
mic variance-limited measurement of the large-scale E-mode 
power spectrum will improve on the current WMAP and 
Planck measurements of the optical depth r to reionization 
by a factor of four. Not only can optical depth measurements 
constrain cosmic reionization scenarios, but it can also further 
constrain the amplitude of primordial density fluctuations A s , 
whose precision is limited by its degeneracy with r. With 
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a tightly constrained A s measurement, comparisons with the 
amplitude of large scale structure of the Universe, erg, can be 
used to infer the sum of neutrino masses to 0.015 eV (Allison 
et al. 2015). Finally, the observed anomalies in the temper¬ 
ature fluctuations on large angular scales require polarization 
measurements on the same scales to probe physics beyond 
standard ACDM (Bennett et al. 2011; Planck Collaboration 
2015d). 
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APPENDIX 

A. SIGNAL-TO-NOISE ON LARGE ANGULAR SCALES 

The purpose of this Appendix is to demonstrate that in the 
large-noise limit, the signal-to-noise of B-mode polarization 
measurements increases with decreasing multipole l. Cosmic 
variance, the uncertainty in the observed CMB anisotropy due 
to it being a single realization of a Gaussian random field, 
is a limitation for large scale experiments in particular since 
the uncertainty in the power spectrum scales as 7“ 1/2 in the 
cosmic variance limit. For a beam window function bi and 
noise power spectrum Nf, the signal-to-noise per multipole of 
a power spectrum is (Knox 1995) 

Q b] = V/ sky (l+l/2) 

For C^b 2 » Nt, i.e. in the cosmic variance limit, the 
l 112 factor dominates the signal-to-noise, but when the power 
spectra of the noise and CMB are comparable the relationship 
is not as simple, and the projected significance can sometimes 
increase at large angular scales. 

CMB polarization experiments are far from the cosmic vari¬ 
ance measurement limits of the B-mode power spectrum, so 
at large angular scales the signal-to-noise in Equation Al be¬ 
comes 


C^V/^+l/2) 

N ( 


(A2) 


consistent with noise, and does not indicate any intrinsic vari¬ 
ation. To estimate the spectral index of thermal dust emission, 
the Planck team subdivided mid-latitude data into 400 over¬ 
lapping patches of sky (Planck Int XXII). They subtracted 
lower frequency data from dust-dominated bands to remove 
the achromatic CMB component and then took ratios of the 
residual at two different frequencies, v\ and V 2 , denoted 

n / x Iv 2 — /vo 

Ry 0 (V2,Vl) = - ~ - 

l v x L 0 

^ lv2_\^ B V1 (T D ) 

\vj B V1 (T D ) 

where fiv is the spectral index of the dust, 7 d is the dust tem¬ 
perature, and I y is the total signal intensity at frequency v. 
To remove the CMB contribution from the dust-dominated 
bands, they subtracted I VQ from I V2 and 7 Vl before taking 
a ratio, and assumed the dust follows a modified black- 
body spectrum oc v^ d B v (T d ). They estimated the dust tem¬ 
perature using R t ( 3000, 857) with DIRBE data as the high- 
frequency template, while they computed spectral indices us¬ 
ing 7 ?j T 0 q } ( 353, 217). Using 7 d, they implicitly solved for y# D 
for each patch. 

The resulting distribution of / 3 d from Planck has a l<x dis¬ 
persion of crg D = 0.17, which at face value suggests true on- 
sky spectral index variation. However, each sky patch has 
an associated local dispersion cr 353 in the 353 GHz dust tem¬ 
plate brightness, with higher values corresponding to higher 
signal-to-noise, and Planck Int XXII claim that low signal-to- 
noise values with cr 353 < 20 yuK are dominated by instrument 
noise. Additionally, in Appendix B of Planck Int XXII, Monte 
Carlo simulations of dust polarization with constant T D and 
J3j) and instrument noise give a ler dispersion crg D = 0.07 in 
J3j) (which is therefore entirely due to instrument noise) for 
sky patches with signal <Hj 53 > 30 /i K. 

Using data from the top plot of Figure 8 of Planck Int XXII 
and the relation from Equation B2 fixing Td to be constant, we 
rescale R^ 00 (353,217), assuming a log-linear relation between 
R^ 00 (353,217) and yS D - We rescale to have the quoted 
mean and standard deviation of (3v = 1.59 + 0.17. We re¬ 
produce Figure 8 of Planck Int XXII in our Figure 7. Us¬ 
ing these estimates, we compute the standard deviation of 
subsets of the data with cr v 353 > [20, 30,40] yuK, obtaining 
<t^ d = [0.12,0.07,0.05], respectively. We observe that the 
observed variation of <x^ d = 0.07 for cr!^ > 30 yuK is entirely 
consistent with effects of instrument noise seen in the Planck 
Int XXII Monte Carlo simulations and, therefore, with the hy¬ 
pothesis of a constant dust spectral index in the polarization 
data. Based on this, we claim that our model of dust spatial 
spectral variation with <x^ d = 0.01 as shown in Figure 3 is 
consistent with the latest Planck data. 
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